#pragma once
#include "fluid.h"
#include "pipe.h"
#include <vector>

class MultiphaseFlowSolver {
public:
    struct Results {
        std::vector<double> pressure;     // Pa
        std::vector<double> temperature;  // K
        std::vector<double> gasHoldup;    // -
        std::vector<double> oilHoldup;    // -
        std::vector<double> waterHoldup;  // -
        std::vector<double> gasVelocity;  // m/s
        std::vector<double> oilVelocity;  // m/s
        std::vector<double> waterVelocity;// m/s
    };

    MultiphaseFlowSolver(const Fluid& fluid, const PipeParameters& params);
    Results solve();

private:
    // 基本参数
    const Fluid& fluid_;
    PipeParameters params_;
    int numGrids_;
    double dx_;
    double dt_;
    
    // 相场变量
    std::vector<double> alpha_g_;     // 气相体积分数
    std::vector<double> alpha_d_;     // 液滴体积分数
    std::vector<double> alpha_l_;     // 连续液相体积分数
    
    // 速度场
    std::vector<double> V_g_;         // 气相速度
    std::vector<double> V_d_;         // 液滴速度
    std::vector<double> V_l_;         // 连续液相速度
    
    // 压力和温度
    std::vector<double> pressure_;    // 压力场
    std::vector<double> pressure_correction_; // 压力修正量
    std::vector<double> temperature_; // 温度场
    
    // 物性参数
    std::vector<double> gasDensity_;    // 气相密度
    std::vector<double> oilDensity_;    // 油相密度
    std::vector<double> waterDensity_;  // 水相密度
    std::vector<double> gasViscosity_;  // 气相粘度
    std::vector<double> oilViscosity_;  // 油相粘度
    std::vector<double> waterViscosity_;// 水相粘度
    
    // 通量存储
    std::vector<double> gas_flux_old_;   // 上一时间步气相通量
    std::vector<double> liquid_flux_old_;// 上一时间步液相通量
    std::vector<double> energy_flux_old_;  // 上一时间步的能量通量
    
    // 上一时间步的值
    std::vector<double> alpha_g_old_;
    std::vector<double> alpha_d_old_;
    std::vector<double> V_g_old_;
    std::vector<double> V_d_old_;
    
    // 主要求解函数
    void solveTimeStep();
    void solveGasDropletMomentum();
    void solvePressureCorrection();
    void solveMassConservation();
    void solveEnergyEquation();
    
    // 初始化和更新函数
    void initializeFlow();
    void initializeArrays(Results& results);
    void updateProperties();
    void updateVelocities(int i);
    void updateVelocityWithPressureCorrection(const std::vector<double>& p_correction);
    
    // 边界条件相关
    void applyBoundaryConditions();
    void applyMassBoundaryConditions();
    void applyTemperatureBoundaryConditions();
    
    // 计算辅助函数
    double calculatePressureGradient(int i);
    void calculatePhasicVelocities(int i);
    void calculateHeatTransfer();
    double calculateFrictionFactor(int i, const std::string& phase);
    double calculateWallShearStress(int i, const std::string& phase);
    double calculateGasWallShear(int i, double Ls, double Lf, double Lu);
    double calculateLiquidWallShear(int i, double Ls, double Lf, double Lu);
    double calculateInterfacialShearStress(int i);
    double calculateSlugLength(double D);
    double calculateUnitLength(double Ls, double Vgs, double Vlf, double Hls, double Hf);
    double calculateFilmLength(double Lu, double Ls);
    double calculateMassImbalance(int i);
    double calculateVelocityCorrectionCoeff(int j);
    
    // 相变和质量传递相关
    double calculateEntrainmentRate(int i);
    double calculateDepositionRate(int i);
    double calculateMassFlowRate(int i);
    double calculateGasMassFlowRate(int i);
    double calculateEquilibriumMassFraction(double p, double T);
    double calculatePhaseChangeRate(int i);
    
    // 数值方法相关
    void solveTridiagonalMatrix(const std::vector<double>& a,
                               const std::vector<double>& b,
                               const std::vector<double>& c,
                               const std::vector<double>& d,
                               std::vector<double>& x);
    
    // 结果处理
    void saveResults(Results& results, double time);
    bool checkConvergence();
};
